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Abstract 

Recently, Chertkov and Chernyak derived an exact expression for the partition sum (nor- 
malization constant) corresponding to a graphical model, which is an expansion around 
the belief propagation (BP) solution. By adding correction terms to the BP free energy, 
one for each "generalized loop" in the factor graph, the exact partition sum is obtained. 
However, the usually enormous number of generalized loops generally prohibits summation 
over all correction terms. In this article we introduce truncated loop series BP (TLSBP), 
a particular way of truncating the loop series of Chertkov & Chernyak by considering gen- 
eralized loops as compositions of simple loops. We analyze the performance of TLSBP in 
different scenarios, including the Ising model on square grids and regular random graphs, 
and on Promedas, a large probabilistic medical diagnostic system. We show that TLSBP 
often improves upon the accuracy of the BP solution, at the expense of increased compu- 
tation time. We also show that the performance of TLSBP strongly depends on the degree 
of interaction between the variables. For weak interactions, truncating the series leads to 
significant improvements, whereas for strong interactions it can be ineffective, even if a 
high number of terms is considered. 

Keywords: Belief propagation, loop calculus, approximate inference, partition function, 
Ising grid, random regular graphs, medical diagnosis 



1. Introduction 

Belief propagation (Pearl, 1988; Murphy et al., 1999) is a popular inference method that 
yields exact marginal probabilities on graphs without loops and can yield surprisingly accu- 
rate results on graphs with loops. BP has been shown to outperform other methods in rather 
diverse and competitive application areas, such as error correcting codes (Gallager, 1963; 
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McEliece et al., 1998), low level vision (Freeman et al., 2000), combinatorial optimization 
(Mezard et al., 2002) and stereo vision (Sun et al., 2005). 

Associated to a probabilistic model is the partition sum, or normalization constant, 
from which marginal probabilities can be obtained. Exact calculation of the partition 
function is only feasible for small problems, and there is considerable statistical physics 
literature devoted to the approximation of this quantity. Existing methods include stochas- 
tic Monte Carlo techniques (Potamianos and Goutsias, 1997) or deterministic algorithms 
which provide lower bounds (Jordan et al., 1999; Leisink and Kappen, 2001), upper bounds 
(Wainwright et al., 2005), or approximations (Yedidia et al., 2005). 

Recently, Chertkov and Chernyak have presented a loop series expansion formula that 
computes correction terms to the belief propagation approximation of the partition sum. 
The series consists of a sum over all so-called generalized loops in the graph. When all 
loops are taken into account, Chertkov & Chernyak show that the exact result is recovered. 
Since the number of generalized loops in a graphical model easily exceeds the number of 
configurations of the model, one could argue that the method is of little practical value. 
However, if one could truncate the expansion in some principled way, the method could 
provide an efficient improvement to BP. 1 

Most inference algorithms on loopy graphs can be viewed as generalizations of BP, 
where messages are propagated between regions of variables. For instance, the junction- 
tree algorithm (Lauritzen and Spiegelhalter, 1988) which transforms the original graph in a 
region tree such that the influence of all loops in the original graph is implicitly captured, 
and the exact result is obtained. However, the complexity of this algorithm is exponential in 
time and space on the size of the largest clique of the resulting join tree, or equivalently, on 
the tree-width of the original graph, a parameter which measures the network complexity. 
Therefore, for graphs with high tree-width one is resorted to approximate methods such 
as Monte Carlo sampling or generalized belief propagation (GBP) (Yedidia et al., 2005), 
which captures the influence of short loops using regions which contain them. One way to 
select valid regions is the cluster variation method (CVM) (Pelizzola, 2005). In general, 
selecting a good set of regions is not an easy task, as described in (Welling et al., 2005). 
Alternatively, double-loop methods can be used (Heskes et al., 2003; Yuille, 2002) which 
are guaranteed to converge, often at the cost of more computation time. 

In this work we propose TLSBP, an algorithm to compute generalized loops in a graph 
which are then used for the approximate computation of the partition sum and the single- 
node marginals. The proposed algorithm is parametrized by two arguments which are used 
to prune the search for generalized loops. For large enough values of these parameters, all 
generalized loops present in a graph are retrieved and the exact result is obtained. One 
can then study how the error is progressively corrected as more terms are considered in the 
series. For cases were exhaustive computation of all loops is not feasible, the search can 
be pruned, and the result is a truncated approximation of the exact solution. We focus 
mainly on problems where BP converges easily, without the need of damping or double 
loop alternatives (Heskes et al., 2003; Yuille, 2002) to force convergence. It is known that 
accuracy of the BP solution and convergence rate are negatively correlated. Throughout 
the paper we show evidence that for those cases where BP has difficulties to converge, loop 

1. Note that the number of generalized loops in a finite graph is finite, and strictly speaking, the term 
series denotes an infinite sequence of terms. For clarity, we prefer to use the original terminology. 
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corrections are of little use, since loops of all lengths tend to have contributions of similar 
order of magnitude. 

The paper is organized as follows. In Section 2 we briefly summarize the series expansion 
method of Chertkov and Chernyak. In Section 3 we provide a formal characterization of 
the different types of generalized loops that can be present in an arbitrary graph. This 
description is relevant to understand the proposed algorithm described in Section 4. We 
present experimental results in Section 5 for the Ising model on grids, regular random graphs 
and medical diagnosis. Concering grids and regular graphs, we show that the success of 
restricting the loop series expansion to a reduced quantity of loops depends on the type 
of interactions between the variables in the network. For weak interactions, the largest 
correction terms come from the small elementary loops and therefore truncation of the 
series at some maximal loop length can be effective. For strong interactions, loops of all 
lengths contribute significantly and truncation is of limited use. We numerically show that 
when more loops are taken into account, the error of the partition sum decreases and when 
all loops are taken into account the method is correct up to machine precision. We also apply 
the truncated loop expansion to a large probabilistic medical diagnostic decision support 
system (Wiegerinck et al., 1999). The network has 2000 diagnoses and about 1000 findings 
and is intractable for computation. However, for each patient case unobserved findings and 
irrelevant diagnoses can be pruned from the network. This leaves a much smaller network 
that may or may not be tractable depending on the set of clamped findings. For a number 
of patient cases, we compare the BP approximation and the truncated loop correction. We 
show results and characterize when the loop corrections significantly improve the accuracy 
of the BP solution. Finally, in Section 6 we provide some concluding remarks. 

2. BP and the Loop Series Expansion 

Consider a probability model on a set of binary variables Xi = ±1, i = 1, . . . , n: 

1 m m 

P(x) = -l[f a (x a ), Z = J2HMXa), (1) 

a=l x a=l 

where a = 1, . . . ,m labels interactions (factors) on subsets of variables x a , and Z is the 
partition function, which sums over all possible states or variable configurations. Note that 
the only restriction here is that variables are binary, since arbitrary factor nodes are allowed, 
as in (Chertkov and Chernyak). 

The probability distribution in (1) can be directly expressed by means of a factor graph 
(Kschischang et al., 2001), a bipartite graph where variable nodes % are connected to factor 
nodes a if and only if x\ is an argument of f a . Figure 3 (left) on page 14 shows an example 
of a graph where variable and factor nodes are indicated by circles and squares respectively. 

For completeness, we now briefly summarize Pearl's belief propagation (BP) (Pearl, 
1988) and define the Bethe free energy. If the graph is acyclic, BP iterates the following 
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message update equations, until a fixed point is reached: 
variable i to factor a: /ij_> Q (xj) 

factor a to variable i: 

where i € a denotes variables included in factor a, and a 3 i denotes factor indices a which 
have i as argument. After the fixed point is reached, exact marginals and correlations 
associated with the factors ("beliefs") can be computed using: 

bi(xi) oc JJ/x a ^j(xi), 

where oc indicates normalization so that beliefs sum to one. 

For graphs with cycles the same update equations can be iterated (the algorithm is 
then called loopy, or iterative, belief propagation), and one can still obtain very accurate 
approximations of the beliefs. However, convergence is not guaranteed in these cases. For 
example, BP can get stuck in limit cycles. An important step towards the understanding 
and characterization of the convergence properties of BP came from the observation that 
fixed points of this algorithm correspond to stationary points of a particular function of the 
beliefs, known as the Bethe free energy (Yedidia et al., 2001), which is defined as: 

Fbp = Ubp — Hbp, (2) 
where Ubp is the Bethe average energy: 

m 

U B p = ~^2^2 b a {x a ) log f a (x a ), 

a=l X a 

and Hpp is the Bethe approximate entropy: 

m n 

H B p = - ^j^jb a (x a )\ogb a (x a ) - 1) ^ bi(xi) log6j(xj). (3) 

Q = l X a 1 = 1 Xi 

where di is the number of neighboring factor nodes of variable node i. The second term 
in (3) ensures that every node in the graph is counted once, see (Yedidia et al., 2005) for 
details. The BP algorithm tries to minimize (2) and, for trees, the exact partition function 
can be obtained after the fixed point has been reached, Z = exp(—Fpp)- However, for 
graphs with loops Fpp provides just an approximation. 

If one can calculate the exact partition function Z defined in Equation (1), one can also 
calculate any marginal in the network. For instance, the marginal 

m 

where Z(9 t ) := £ e e ^ J] f a (x a ) 

8i^0 x a=l 



= n vp^iixi), 

I35i\{a} 

x a \{i} jea\{i} 



Pi(Xi) 



dlog Z(0i 



d6i(xi) 
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is the partition sum of the network, perturbed by an additional local field potential 9i on 
variable X{. 

Alternatively, one can compute different partition functions for different settings of the 
variables, and derive the marginals from ratios of them: 

where Z Xi indicates the partition function calculated from the same model conditioning 
on variable i, i.e, with variable i fixed (clamped) to value X{. Therefore, approximation 
errors in the computation of any marginal can be related to approximation errors in the 
computation of Z. We will thus focus on the approximation of Z mainly, although marginal 
probabilities will be computed as well. 

Of central interest in this work is the concept of generalized loop, which is defined in 
the following way: 

Definition 1 A generalized loop in a graph G = (V,E) is any subgraph C = (V',E r ), 
V ^V,E'C_ (V x V') n E such that each node in V' has degree two or larger. The length 
(size) of a generalized loop is its number of edges. 

For the rest of the paper, the terms loop and generalized loop are used interchangeably. 
The main result of (Chertkov and Chernyak) is the following. Let b a (x a ), h(xi) denote the 
beliefs after the BP algorithm has been converged, and let Z^p = exp(— Fbp) denote the 
corresponding approximation to the partition sum, with Fbp the value of the Bethe free 
energy evaluated at the BP solution. Then Zp,p is related to the exact partition sum Z as: 



Z = Z BP ( 1 + £ r (C)\ , r{C) = H ^(C) [J M«(C) 

V cec / i€C aec 



(5) 



where summation is over the set C of all generalized loops in the factor graph. Any term 
r(C) in the series corresponds to a product with as many factors as nodes present in the 
loop. Each factor is related to the beliefs at each variable node or factor node according to 
the following formulas: 

, (r) (1-tth)^- 1 + {-1)^(1 + rn i )^- 1 n(r) V 1 ffil 

^ C >- 2 ( 1 - m 2 ^(C)-i ' U^)- 2^ 1 W 

Ma(C) = /^ J b a (x a ) [j (xi-rrii), (7) 

x a ieC,iea 

where nii = Yl x = ^«(+) — ) is the expected value of Xi computed in the BP 

approximation. Generally, terms r{C) can take positive or negative values. Even the same 
variable i may have positive or negative subterms fii depending on the structure of the 
particular loop. 

Expression (5) represents an exact and finite decomposition of the partition function 
with the first term of the series being exactly represented by the BP solution. Note that, 
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although the series is finite, the number of generalized loops in the factor graph can be 
enormous and easily exceed the number of configurations 2 n . In these cases the loop series 
is less efficient than the most naive way to compute Z exactly, namely by summing the 
contributions of all 2™ configurations one by one. 

On the other hand, it may be that restricting the sum in (5) to a subset of the total 
generalized loops captures the most important corrections and may yield a significant im- 
provement in comparison to the BP estimate. We therefore define the truncated form of 
the loop corrected partition function as: 



Ztlsbp = Z BP 1 1 + ^ r(C) J 



(8) 



where summation is over the subset C C C obtained by Algorithm 2, which we will discuss 
in Section 4. Approximations for the single- node marginals can then be obtained from (8), 
using the method proposed in Equation (4): 

/ j ^TLSBP 

Because the the terms r{C) can have different signs, the approximation Ztlsbp is in general 
not a bound of the exact Z, but just an approximation. 



3. Loop Characterization 

In this section we characterize different types of generalized loops that can be present in 
a graph. This classification is the basis of the algorithm described in the next section and 
also exemplifies the different shapes a generalized loop can take. For clarity, we illustrate 
them by means of a factor graph arranged in a square lattice with only pairwise interactions. 
However, definitions are not restricted to this particular model and can be applied generally 
to any factor graph. 

Definition 2 A simple (elementary) generalized loop (from now on simple loop) is defined 
as a connected subgraph of the original graph where all nodes have exactly degree two. 

This type of generalized loop coincides with the concept of simple circuit or simple cycle 
in graph theory: a path which starts and ends at the same node with no repeated vertices 
except for the start and end vertex. Figure la shows an example of a simple loop of size 8. 
On the contrary, in Figure lb we show an example of generalized loop which is not a simple 
loop, because three nodes have degree larger than two. 

We now define the union of two generalized loops, l\ = (V±,Ei) and li = (V^-E^), as 
the generalized loop which results from taking the union of the vertices and the edges of l\ 
and I2, that is, I' = l\ U I2 = (V\ U V2, E\ U E2). Note that the union of two simple loops 
is never a simple loop except for the trivial case in which both loops are equal. Figure lb 
shows an example of a generalized loop which can be described as the union of three simple 
loops, each of size 8. The same example can be also defined as the union of two overlapping 
simple loops, each of size 12. 
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Figure 1: Examples of generalized loops in a factor graph with lattice structure. Variable 
nodes and factor nodes are represented as squares and rhombus respectively. 
Generalized loops are indicated using bold edges underlying the factor graph, 
(a) A simple loop, (b) A non-simple loop which is neither a disconnected loop 
nor a complex loop, (c) A disconnected loop of three components, each a simple 
loop, (d) A disconnected loop of two components, the left one a non-simple loop, 
(e) A complex loop which is not a disconnected loop, (f) A complex loop which 
is also a disconnected loop. (See text for definitions). 
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Figure 2: Diagrammatic representation of the different types of generalized loops present 
in any graph. Sizes of the sets are just indicative and depend on the particular 
instance. 



Definition 3 A disconnected generalized loop, disconnected loop, is defined as a gener- 
alized loop with more than one connected component. 

Figure lc shows an example of a disconnected loop composed of three simple loops. 
Note that components are not restricted to be simple loops. Figure Id illustrates this fact 
using an example where one connected component (the left one) is not a simple loop. 

Definition 4 A complex generalized loop, complex loop, is defined as a generalized loop 
which cannot be expressed as the union of two or more different simple loops. 

Figures le and If are examples of complex loops. Intuitively, they result after the 
connection of two or more connected components of a disconnected loop. 

Any generalized loop can be categorized according to these three different categories: 
a simple loop cannot be a disconnected loop, neither a complex loop. On the other hand, 
since Definitions 3 and 4 are not mutually exclusive, a disconnected loop can be a complex 
loop and vice- versa, and also there are generalized loops which are neither disconnected nor 
complex, for instance the example of Figure lb. An example of a disconnected loop which 
is not a complex loop is shown in Figure lc. An example of a complex loop which is not 
a disconnected loop is shown in Figure le. Finally, an example of a complex loop which is 
also a disconnected loop is shown in Figure If. 

We finish this characterization using a diagrammatic representation in Figure 2 which 
illustrates the definitions. Usually, the smallest subset contains the simple loops and both 
disconnected loops and complex loops have nonempty intersection. There is another subset 
of all generalized loops which are neither simple, disconnected, nor complex. 
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4. The Truncated Loop Series Algorithm 

In this section we describe the TLSBP algorithm to compute generalized loops in a factor 
graph, and use them to correct the BP solution. The algorithm is based on the principle 
that every generalized loop can be decomposed in smaller loops. The general idea is to 
search first for a subset of the simple loops and, after that, merge all of them iteratively 
until no new loops are produced. As expected, a brute force search algorithm will only 
work for small instances. We therefore prune the search using two different bounds as 
input arguments. Eventually, a high number of generalized loops which presumably will 
account for the major contributions in the loops series expansion will be obtained. We show 
that the algorithm is complete, or equivalently, that all generalized loops are obtained by 
the proposed approach when the constraints expressed by the two arguments are relaxed. 
Although exhaustive enumeration is of little interest for complex instances, it allows to 
check the validity of (5) and to study the loop series expansion for simpler instances. The 
algorithm is composed of three steps: 

1. First, we remove recursively all the leaves of the original graph, until its 2-core is 
obtained. This initial step has two main advantages. On the one hand, since some 
nodes are deleted, the complexity of the problem is reduced. On the other hand, we 
can use the resulting graph as a test for any possible improvement to the BP solution. 
Indeed, if the original graph did not contain any loop then the null graph is obtained, 
the BP solution is exact on the original graph, and the series expansion has only one 
term. On the other hand, if a nonempty graph remains after this preprocessing, it 
will have loops and the BP solution can be improved using the proposed approach. 

2. After the graph is preprocessed, the second step searches for simple loops. The result of 
this search will be the initial set of loops for the next step and will also provide a bound 
b which will be used to truncate the search for new generalized loops. Finding circuits 
in a graph is a problem addressed for long (Tiernan, 1970; Tarjan, 1973; Johnson, 
1975) whose computational complexity grows exponentially with the length of the 
cycle (Johnson, 1975). Nevertheless, we do not count all the simple loops but only 
a subset. Actually, to avoid dependence on particular instances, we parametrize this 
search by a size S, which limits the number of shortest simple loops to be considered. 
Once S simple loops have been found in order of increasing length, the length of the 
largest simple loop is used as the bound b for the remaining steps. 

3. The third step of the algorithm consists of obtaining all non-simple loops that the set 
of S simple loops can "generate" . 

According to definition 4, complex loops can not be expressed as union of simple loops. 
To develop a complete method, in the sense that all existing loops can be obtained, we 
define the operation merge loops, which extends the simple union in such a way that complex 
loops are retrieved as well. Given two generalized loops, l\,l2, merge loops returns a set 
of generalized loops. One can observe that for each disconnected loop, a set of complex 
loops can be generated by connecting two (or more) components of the disconnected loop. 
In other words, complex loops can be expressed as the union of disjoint loops with a path 
connecting two vertices of different components. Therefore the set computed by merge loops 



9 



ViCENg Gomez and J.M. Mooij and H.J. Kappen 



will have only one element I' = {l\ VJl 2 } if Zi U l 2 is not disconnected. Otherwise, all the 
possible complex loops in which l\ U l 2 appears are included in the resulting set. 

We use the following procedure to compute all complex loops associated to the discon- 
nected loop V: we start at a vertex of a connected component of V and perform depth-first- 
search (DFS) until a vertex of a different component has been reached. At this point, the 
connecting path and the reached component are added to the first component. Now the 
generalized loop has one less connected component. This procedure is repeated again until 
the resulting generalized loop is not disconnected, or equivalently, until all its vertices are 
members of the first connected component. Iterating this search for each vertex every time 
two components are connected, and also for each initial connected component, one obtains 
all the required complex loops. 

Algorithm 1 merge loops 
Arguments: 

h = (Vi,£i) loop, 

h = (V 2 ,E 2 ) loop, 

b maximal length of a loop, 

M maximal depth of complex loops search, 

G preprocessed factor graph 

1: newloops <— 

2: if (\E 1 U E 2 \ < b) then 

3: C <— Find connected components^! U l 2 ) 

4: newloops <— {l\ U l 2 } 

5: for all (a e C) do 

6: for all (vi <G q) do 

7: newloops <— newloops U Find complex loopsDFS(i>j, q, C, M, b, G) 

8: end for 
9: end for 
10: end if 

11: return newloops 



Note that deciding whether l\ U l 2 is disconnected or not requires finding all connected 
components of the resulting loop. Moreover, given a disconnected loop, the number of 
associated complex loops can be enormous. In practice, the bound b obtained previously 
is used to reduce the number of calculations. First, testing if the length of l\ U l 2 is larger 
than b can be done without computing the connected components. Second, the DFS search 
for complex loops is limited using b, so very large complex loops will not be retrieved. 

However, restricting the DFS search for complex loops using the bound b could result 
in too deep searches. Consider the worst case of merging the two shortest, non-overlapping, 
simple loops which have size L s . The maximum depth of the DFS search for complex loops 
is d = b — 2L S . Then the computational complexity of the merge loops operation depends 
exponentially on d. This dependence is especially relevant when b » L s , for instance 
in cases where loops of many different lengths exist. To overcome this problem we define 
another parameter M, the maximum depth of the DFS search in the merge loops operation. 
For small values of M, the operation merge loops will be fast but a few (if any) complex 
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loops will be obtained. Conversely, for higher values of M the operation merge loops will 
find more complex loops at the cost of more time. 

Algorithm 1 in the previous page describes briefly the operation merge loops. It receives 
two loops h and I2, and bounds b and M as arguments, and returns the set newloops which 
contains the loop resulting of the union of l\ and I2 plus all complex loops obtained in the 
DFS search bounded by b and M. 

Algorithm 2 Algorithm TLSBP 
Arguments: 

S maximal number of simple loops, 

M maximal depth of complex loops search, 

G original factor graph 

1: Run belief propagation algorithm over G 
2: G' <- Obtain the 2-core(G) 
3: C <- 

4: if (-iempty(G')) then 

5: (sloops, b) <— Compute first S simple loops(G') 
6: {oldloops, newloops) <— (sloops,®) 
7: C <— sloops 

8: while (-iempty (oldloops)) do 
9: for all (Zi € sloops) do 
10: for all (I2 £ oldloops) do 

11: newloops <— newloops U mergeLoops(/i, I2, b, M, G') 

12: end for 

13: end for 

14: oldloops <— newloops 

15: C f- C U newloops 

16: end while 

17: end if 

18: return the result of expression (8) using C 



Once the problem of expressing all generalized loops as compositions of simple loops 
has been solved using the merge loops operation, we need to define an efficient procedure to 
merge them. Note that, given S simple loops, a brute force approach tries all combinations 
of two, three, ... up to S 1 — 1 simple loops. Hence the total number is: 

which is prohibitive. Nevertheless, we can avoid redundant combinations by merging pairs 
of loops iteratively: in a first iteration, all pairs of simple loops are merged, which produces 
new generalized loops. In a next iteration i, instead of performing all J mergings, only the 
new generalized loops obtained in iteration i — 1 are merged with the initial set of simple 
loops. The process ends when no new loops are found. Using this merging procedure, 
although the asymptotic cost is still exponential in S, many redundant mergings are not 
considered. 
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Summarizing, the third step applies iteratively the merge loops operation until no new 
generalized loops are obtained. After this step has finished, the final step computes the 
truncated loop corrected partition function defined in Equation (8) using all the obtained 
generalized loops. We describe the full procedure in Algorithm 2. Lines 2 and 4 correspond 
to the first and second steps and lines 5 — 13 correspond to the third step. 

To show that this process produces all the generalized loops we first assume that S is 
sufficiently large to account for all the simple loops in the graph, and that M is larger or 
equal than the number of edges of the graph. Now let C be a generalized loop. According 
to the definitions of Section 3, either C can be expressed as a union of s simple loops, or C 
is a complex loop. In the first case, C is clearly produced in the sth iteration. In the second 
case, let s' denote the number of simple loops which appear in C. Then C is produced in 
iteration s', during the DFS for complex loops within the merging of one of the simple loops 
contained in C. 

The obtained collection of loops can be used for the approximation of the singe node 
marginals as well, as described in Equation (9). The method consists of clamping one 
variable i to all its possible values (±1) and computing the corresponding approximations 
of the partition functions: Ztlsbp and Z^Ysbp- This requires to run BP in each clamped 
network, and reuse the set of loops replacing with zero those terms where the clamped 
variable appears. The computational complexity of approximating all marginals using this 
approach is in general 0(N ■ L • d ■ Tbp), where L is the number of found loops, d is the 
cardinality of the variables (two in our case), and Tbp the average time of BP to converge 
after clamping one variable. Usually, this task requires less computation time than the 
search for loops. 

As a final remark, we want to stress a more technical aspect related to the implemen- 
tation. Note that generalized loops can be expressed as the composition of other loops in 
many different ways. In consequence, they all must be stored incrementally and the opera- 
tion of checking if a loop has been previously counted or not should be done efficiently. An 
appropriate way to implement this fast look- up/insertion is to encode all loops in a string 
composed by the edge identifiers in some order with a separator character between them. 
This identifier is used as a key to index an ordered tree, or hash structure. In practice, a 
hash structure is only necessary if large amounts of loops need to be stored. For the cases 
analyzed here, choosing a balanced tree instead of a hash table resulted in a more efficient 
data structure. 

5. Experiments 

In this section we show the performance of TLSBP in two different scenarios. First, we focus 
on square lattices and study how loop corrections improve the BP solution as a function 
of the interaction between variables and the size of the problem. Second, we study the 
performance of the method in random regular graphs as a function of the degree between 
the nodes. 

In all the experiments we show results for tractable instances, where the exact solution 
using the junction tree (Lauritzen and Spiegelhalter, 1988) can be computed. Performance is 
evaluated comparing the TLSBP error against the BP solution, and also against the cluster 
variation method (CVM). Instead of using a generalized belief propagation algorithm (GBP) 
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which usually requires several trials to find the proper damping factor to converge, we use 
a double-loop implementation which has convergence guarantees (Heskes et al., 2003). For 
this study we select as outer regions of the CVM method all maximal factors together with 
all loops that consist up to four different variables. This choice represents a good trade-off 
between computation time required for convergence and accuracy of the solution. 

We report two different error measures. Concerning the partition function Z we com- 
pute: 



Error^/ = 



logZ' 



logZ 



(10) 



where Z' is the partition function corresponding to the method used: BP, TLSBP, or 
CVM. Error of single-node marginals was measured using the maximum error, which is 
a reasonable quantity if one is interested in worst-case scenarios: 

Error 6 = max \Pi(xi) - &j(xj)|, (11) 
i=l,...,n 

Xi = ±l 

were again 6j(xj) are the single-node marginal approximations corresponding to the method 
used. 

We use four different schemas for belief- updating of BP: (i) fixed and (ii) random sequen- 
tial updates, (Hi) parallel (or synchronous) updates, and (iv) residual belief propagation 
(RBP), a recent method proposed in (Elidan et al., 2006). The latter method schedules 
the updates of the BP messages heuristically by selecting the next message to be updated 
which has maximum residual, a quantity defined as an upper bound on the distance of the 
current messages from the fixed point. In general, we experienced that for some instances 
where the RBP method converged, the other update schemas (fixed, random sequential and 
parallel updates) failed to converge. 

In all schemas we interpret that a fixed point is reached at iteration t when the maximum 
absolute value of the updates of all the messages from iteration t — 1 to t is smaller than 
a threshold We notice a large correlation between the order of magnitude of $ and the 
ratio between the BP and the TLSBP errors. For this reason we used a very small value of 
the threshold, ■& = 10~ 17 . 



5.1 Ising Grids 

This model is defined on a grid where each variable, also called spin, takes binary values 
Xi = ±1. A spin is coupled with its direct neighbors only, so that pairwise interactions 
fij(xi,Xj) = exp(6ijXiXj) are considered, parametrized by 9ij. Every spin can be exposed 
to an external field fi(xi) = exp(6>jXj), or single-node potential, parametrized by Q{. Figure 3 
(left) shows the factor graph associated to the 4x4 Ising grid, composed of 16 variables. The 
Ising grid model is often used as a test-bed for inference algorithms. It is of great relevance 
in statistical physics, and has applications in different areas such as image processing. In 
our context it also represents a challenge since it has many loops. Good results in this 
model will likely translate into good results for less loopy graphs. 

Usually, two cases are differentiated according to the sign of the 9ij parameters. For 
6ij > coupled spins tend to be in the same state. This is known as the attractive, 
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or "ferromagnetic" setting. On the other hand, for mixed interactions, 9^ can be either 
positive or negative, and this setting is called "spin-glass" configuration. Concerning the 
external field, one can distinguish two cases. For the case of nonzero fields, larger values 
of 9i imply easier inference problems in general. On the other hand, for Oi = 0, there exist 
two phase transitions from easy inference problems (small %) to more difficult ones (large 
Oij) depending on the type of pairwise couplings (see Mooij and Kappen, 2005, for more 
details). 

This experimental subsection is structured in three parts: First, we study a small 4x4 
grid. We then study the performance of the algorithm in a 10x10 grid, where complete 
enumeration of all generalized loops is not feasible. Finally, we analyze the scalability of 
the method with problem size. 

The 4x4 Ising grid is complex enough to account for all types of generalized loops. It is 
the smallest size where complex loops are present. At the same time, the problem is still 
tractable and exhaustive enumeration of all the loops can be done. 

We ran the TLSBP algorithm in this model with arguments S and M large enough to 
retrieve all the loops. Also, the maximum length b was constrained to be 48, the total num- 
ber of edges for this model. After 4 iterations all generalized loops were obtained. The total 
number is 16371 from which 213 are simple loops. The rest of generalized loops are classified 
as follows: 174 complex and disconnected loops, 1646 complex but non-disconnected loops, 
604 non-complex but disconnected loops, and 13734 neither complex nor disconnected loops. 

Figure 3 (right) shows the histogram of all generalized loops for this small grid. Since we 
use the factor graph representation the smallest loop has length 8. The largest generalized 
loop includes all nodes and all edges of the preprocessed graph, and has length 48. The 
Poisson-like shape of the histogram is a characteristic of this model and for larger instances 




Figure 3: (left) A factor graph representing the 4x4 Ising grid, (right) Number of gener- 
alized loops as a function of the length using the factor graph representation. 



14 



Truncating the Loop Series Expansion for BP 



we observed the same tendency. Thus the analysis for this small model can be extrapolated 
to some extent to grids with more variables. 

To analyze how the error changes as more loops are considered it is useful to sort all 
the terms r(C) by their absolute value in descending order such that |r(Cj)| > |r(Cj+i)|. 
We then compute, for each number of loops I = 1... 16371, the approximated partition 
function which accounts for the I most important loops: 



From these values of the partition function we calculate the error measure indicated in 
Equation (11). Estimations of the single-node marginals were obtained using the clamping 
method, and their corresponding error was calculated using Equation (10). 

We now study how loop contributions change as a function of the coupling strength 
between the variables. We ran several experiments using mixed interactions with 6ij ~ 
A/"(0, a 2 ) independently for each factor node, and a varying between 0.1 and 2. Single-node 
potentials were drawn according to 6i ~ Af(0, 0.05 2 ). For small values of a, interactions are 
weak and BP converges easily, whereas for high values of a variables are strongly coupled 
and BP has more difficulties, or does not converge at all. 

Figure 4 shows results of representative instances of three different interaction strengths. 
For each instance we plot the partition function error (left column) together with errors in 
the single-node marginals (middle column) and loop contributions as a function of the length 
(right column). First, we can see that improvements of the partition sum correspond to 
improvements of the estimates of marginal probabilities as well. Second, for weak couplings 
(<t = 0.1, first row) we can see that truncating the series until a small number of loops 
(around 10) is enough to achieve machine precision. In this case the errors in BP are most 
prominently due to small simple loops. As the right column illustrates, loop contributions 
decrease exponentially with the size, and loops with the same length correspond to very 
similar contributions. Larger loops give negligible contributions and can thus be ignored 
by truncating the series. As interactions are strengthened, however, more loops have to 
be considered to achieve maximum accuracy, and contributions show more variability for a 
given length (see middle row). Also, oscillations of the error due to the different signs in 
loop terms (caused by the mixed interactions) of the same order of magnitude become more 
frequent. Eventually, for large couplings (a > 2), where BP often fails to converge, loops 
of all lengths give significant contributions. In the bottom panels of Figure 4 we show an 
example of a 'difficult' case for which the BP result is not improved until more than 10 3 
loop terms are summed. The observed discontinuities in the error of the partition sum are 
caused by the fact that oscillations become more pronounced, and corrections composed 
of negative terms rj(Cj) can result in negative values of the partially corrected partition 
function, see Equation (12). This occurs for very strong interactions only, and when a 
small fraction of the total number of loops is considered. In addition, as the right column 
indicates, there is a shift of the main contributions towards the largest loops. 

After analyzing a small grid, we now address the case of the 10x10 Ising grid, where ex- 
haustive enumeration of all the loops is not computationally feasible. We test the algorithm 
in two scenarios: for attractive interactions (ferromagnetic model) where pairwise interac- 




(12) 



15 



ViCENg Gomez and J.M. Mooij and H.J. Kappen 



Coupling strength o = 0.1 



Coupling strength o = 0.1 



Coupling strength o = 0.1 



10" 



N 10 



m io- 10 



10" 



10" 



N 10 



m io" 10 



- Error BP 

- Error TLSBP 



10° 10 1 10 2 10 3 10 4 



#loops 
Coupling strength a = 0.5 



10" 




10° 10 1 10 2 10 3 10 4 
#loops 

Coupling strength a = 2 



10" 



|10" 5 



2 10 

LU 



- Error BP 

- Error TLSBP 



10 



10" 



10° 10 1 10 2 10 3 10 4 
#loops 

Coupling strength a = 0.5 



■§,10" 5 



2 10" 



10" 




10 



9. 10" 



10" 



10 



10" 



9. io" 



10 







10 20 30 40 50 
Loop length 

Coupling strength a = 0.5 



10° 10 1 10 2 10 3 10 4 
#loops 

Coupling strength o = 2 



N 10 




10" 



10 



10 20 30 40 50 
Loop length 

Coupling strength a = 2 



10° 10 1 10 2 10 3 10 
#loops 



10" 



10" 



#loops 



10 20 30 40 50 
Loop length 



Figure 4: Cumulative error for the spin-glass 4x4 Ising grid for different interaction 
strengths, see Equation (12). (left column) Error of Z. (middle column) 
Error of single-node marginals. Dashed lines correspond to the BP error, and 
solid lines correspond to the loop-corrected (TLSBP) error, (right column) 
Absolute values of all loop terms as a function of the length of the corresponding 
loop. 



tions are parametrized as 9ij = |#^-|, 0'^ ~ J\f(0, a 2 ), and also for the previous case of mixed 
interactions (spin-glass model). Single-node potentials were chosen 9i ~ M(0.l, 0.05 2 ) in 
both cases. 
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Attractive interactions 
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Figure 5: TLSBP error for the 10x10 Ising grid with attractive interactions for different 
values of the parameter S. (left) Error of the partition function, (right) Error 
of single-node marginals. 



We show results in Figures 5 and 6 for three values of the parameter S = {10, 100, 1000} 
and a fixed value of M = 10. For S = 10 and S = 100 only simple loops were obtained 
whereas for S = 1000 a total of 44590 generalized loops was used to compute the truncated 
partition sum. Results are averaged errors over 50 random instances. The selected loops 
were the same in all instances. Although in both types of interactions the BP error (solid 
line with dots) is progressively reduced as more loops are considered, the picture differs 
significantly between the two cases. 

For the ferromagnetic case shown in Figure 5 we noticed that all loops have positive 
contributions, r(C) > 0. This is a consequence of this particular type of interactions, since 
all magnetizations have the same sign at the BP fixed point, and also all nodes have an even 
number of neighbors. Consequently, improvements in the BP result are monotonic as more 
loops are considered, and in this case, the TLSBP can be considered as a lower bound of the 
exact solution. For the case of S = 1000, the BP error is reduced substantially at nonzero a, 
but around a ~ 0.5, where the BP error reaches a maximum, also the TLSBP improvement 
is minimal. From this maximum, the BP error decreases again, and loop corrections tend 
to improve progressively the BP solution again as the coupling is strengthened. We remark 
that improvements were obtained for all instances in the three cases. 

Comparing with CVM, TLSBP is better for weak couplings and for S = 1000 only This 
indicates that for intermediate and strong couplings one would need more than the selected 
44590 generalized loops to improve on the CVM result. 

For the case of spin-glass interactions we report different behavior. From Figure 6 we see 
again that for weak couplings the BP error is corrected substantially, but the improvement 
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Mixed interactions 
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Figure 6: TLSBP error of the 10x10 Ising grid with mixed interactions for different values 
of the parameter S. (left) Error of the partition function, (right) Error of 
single- node marginals. 



decreases as the coupling strength is increased. Eventually, for a ~ 1 BP fails to converge 
in most of the cases and also gives poor results. In these cases loop corrections are of little 
use, and there is no actual difference in considering S = 1000 or S = 10. Moreover, because 
loop terms r(C) now can have different signs, truncating the series can lead to worse results 
for S = 1000 than for S = 10. Interestingly, the range where TLSBP performs better than 
CVM is slightly larger in this type of interactions, TLSBP being better for a < 0.5. 

To end this subsection, we study how loop corrections scale with the number of nodes 
in the graph. We only use spin-glass interactions, since it is a more difficult configuration 
than the ferromagnetic case, as previous experiments suggest. We compare the performance 
for weak couplings [a = 0.1), and strong couplings (cr = 0.5), where BP has difficulties to 
converge in large instances. The number of variables N 2 is increased for grids of size N x N 
until exact computation using the junction tree algorithm is not feasible anymore. 

Since the number of generalized loops grows very fast with the size of the grid, we 
choose increasing values of S as well. We use values of S proportional to the number of 
variable nodes A^ 2 such that S = 10iV 2 . This simple linear increment in S means that as 
N is increased, the proportion of simple loops captured by TLSBP over the total existing 
number of simple loops decreases. It is interesting to see how this affects the performance 
of TLSBP in terms of time complexity and accuracy of the solution. For simplicity, M is 
fixed to zero, so no complex loops are considered. Moreover, to facilitate the computational 
cost comparison, we only compute mergings of pairs of simple loops. Actually, for large 
instances the latter choice does not modify the final set of loops, since generalized loops 
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Figure 7: Scalability of the method in the Ising model, (a) Time complexity as a function 
of the produced number of generalized loops, (b) Relation between the time 
complexity of TLSBP and CVM. Comparison of the error of the partition function 
between BP, TLSBP and CVM as a function of the graph size for (c) weak 
interactions and (d) strong interactions. 



which can only be expressed as compositions of three or more simple loops are pruned using 
the bound b. 

In Figure 7, the top panels show averaged results of the computational cost. The left 
plot indicates the relation between the number of loops computed by TLSBP and the time 
required to compute them. This relation can be fit accurately using a line which means that 
for this choice of parameters S and M, and considering only mergings of simple loops, the 
computational complexity of the algorithm grows just linearly with the found loops. One 
has to keep in mind that the number of loops obtained using the TLSBP algorithm grows 
much faster, but much less than the total number of existing loops in the model. 

Figure 7b shows the CPU time consumed by CVM, TLSBP, and the junction tree 
algorithm. In this case, since we only compute the partition function Z, the CPU time of 
TLSBP is constant for both weak and strong couplings. On the contrary, CVM depends 
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on the type of interactions. For weak interactions, TLSBP is in general more efficient than 
CVM, although the scaling trend is slightly better for CVM. For N = 19, CVM starts to 
be more efficient than TLSBP. For strong interactions, CVM needs significantly more time 
to converge in all cases. If we compare the computational cost of the exact method against 
TLSBP, we can see that the junction tree is very efficient for networks with small N, and the 
best option in those cases. However, for N > 17, the junction tree needs more computation 
time, and for TV > 19, the tree-width of the resulting grids is too large. TLSBP memory 
requirements were considerably less in these cases, since loops can be stored efficiently using 
sets of chars. Also, we can see that the TLSBP scaling is better for this choice of parameters 
than that of the exact method. 

Bottom panels show the accuracy of the TLSBP solution. For weak couplings (bottom- 
left) the BP error is always decreased significantly for this choice of parameters and the 
improvement remains almost constant as N increases, meaning that, in this case the number 
of loops which contribute most to the series expansion does not grow significantly with N. 
Interestingly, results are always better than CVM for this regime. 

For strong couplings (bottom-right) the picture changes. First, results differ more be- 
tween instances causing a less smooth curve. Second, the TLSBP error also increases with 
the problem size, so improvements tend to decrease with N, even faster than the BP error 
decay. Eventually, for the largest tractable instance the TLSBP improvement is still signifi- 
cant, about one order of magnitude. Comparing against CVM, unlike in the weak coupling 
scenario, the TLSBP method does not seem to perform better, and only for some cases 
TLSBP error is comparable to the CVM error on average. The accuracy of the TLSBP 
solution for these instances can be increased by considering larger values of S and M, at 
the cost of more time. 

5.2 Random Graphs 

The previous experimental results were focused on the Ising grid which only considers 
pairwise and singleton interactions in such a way that each node in the graph is at most 
linked with four neighbors. Here we briefly analyze the performance of TLSBP applied on 
a more general case, where interactions are less restricted. 

We perform experiments on random graphs with regular topology, where each variable 
is coupled randomly with d other variables using pairwise interactions parametrized by 
6ij ~ A/"(0, a 2 ). Single-node potentials were parametrized in this case by 8i ~ A/"(0, 0.05 2 ). 
We study how loop corrections improve the BP solution as a function of the degree d, and 
compare improvements against the CVM. As in the previous subsection, for CVM we select 
the loops of four variables and all maximal factors as outer clusters. 

Note that the rate of increase in the number of loops with the degree d is even higher 
than with the number of variables in the Ising model. Adding one more link to all the 
variables means adding N more factor nodes to the factor graph. This raises the number 
of loops dramatically. 

For this scenario, we use A = 20 variables and also increase S every time d is increased. 
We simply start with S = 10 and use increments of 250 for each new d. M was set to 10, and 
all possible mergings were computed. We analyze two scenarios, weak (a = 0.1) and strong 
couplings (a = 0.5), and report averages over 60 random instances for each configuration. 
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As Figure 9 (right) indicates, for a = 0.1 BP converged in all instances, whereas for a = 0.5 
BP convergence becomes more difficult as we increase d. 

Figure 8 (top) shows results for weak interactions. The TLSBP algorithm always corrects 
the BP error, although as d increases, the improvement is progressively reduced. We also 
notice that in all cases and methods the approximation of the partition function (left) 
is less accurate than the approximation of the marginals (right). For d = 15, TLSBP 
improvements are still about one order of magnitude for the partition function, and even 
better for the marginals. As in previous experiments with square lattices, the TLSBP 
approach is generally better than CVM in the weak coupling regime. Here, it is also more 




Figure 8: Results on random regular graphs. TLSBP and CVM errors as a function of the 
degree d. Results are averages over 60 random instances. Errors in the partition 
function for weak interactions (a), marginals for weak interactions (b), partition 
function for strong interactions (c), and marginals for strong interactions (d). 
Insets show percentage of instances where the BP error was corrected. 
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Figure 9: Results on random regular graphs, (a) Computation time of TLSBP and CVM. 

In this case, we averaged also all instances over all 60 weak and 60 strong inter- 
actions, since costs were very similar in both cases, (b) Fraction of the instances 
were BP converged. No convergence is reported when none of the four proposed 
schedules converged. 



stable, since for some dense networks the CVM error can be very large, as we can see for 
d = 13 and d = 15. 

For strong interactions (bottom panels), we see that differences between approximations 
of the partition function and single-node marginals are more remarkable than in the previous 
case. The BP partition function is corrected by TLSBP in more than half of the instances 
for all degrees (see inset of Figure 8c, where we plot the fraction of instances where BP was 
corrected in those cases that converged), although for higher degrees, the TLSBP corrections 
are small using this choice of parameters. On the other hand, single-node BP marginals are 
corrected in almost all cases. In contrast, the CVM approach with our selection of outer 
clusters does not perform better than TLSBP in general. In particular, we see that CVM 
estimates of the partition function are very degraded as networks become more dense. This 
unsatisfactory performance of CVM in the estimation of the partition function is not as 
noticeable in the marginal estimates, where BP results are often improved, although with 
much more variability than the TLSBP method. Interestingly, for those few instances of 
dense networks for which BP converged, CVM estimates of the marginals were very similar 
to TLSBP. 

Finally, we compare computational costs in Figure 9 (left). CVM requires significantly 
more time to converge than the time required by TLSBP searching for loops and calculating 
marginals. If we analyze in detail how the TLSBP cost changes, we can notice different 
types of growth for d < 7 and for d>7. The reason behind these two scaling tendencies can 
be explained by the choice of TLSBP parameters, and the bound b (the size of the largest 
simple loop). For d < 7, many simple loops of different lengths are obtained. Consequently, 
the cost of the merging step grows fast, since many loops with length smaller than b are 
produced. On the other hand, for d > 7 simple loops have similar lengths and, therefore, 
less combinations result in additional loops with length larger than the bound b. Without 
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bounding the length of the loops in the merging step, we would expect the first scaling 
tendency (d < 7) also for values of d > 7. 

From these experiments we can conclude that TLSBP performance is generally better 
than CVM in this domain. We should mention that alternative choices of regions would 
have lead to different CVM results, but will probably not change this conclusion. 

5.3 Medical Diagnosis 

We now study the performance of TLSBP on a "real-world" example, the Promedas medical 
diagnostic network. The diagnostic model in Promedas is based on a Bayesian network. 
The global architecture of this network is similar to QMR-DT (Shwe et al., 1991). It 
consists of a diagnosis layer that is connected to a layer with findings. In addition, there 
is a layer of variables, such as age and gender, that may affect the prior probabilities of 
the diagnoses. Since these variables are always clamped for each patient case, they merely 
change the prior disease probabilities and are irrelevant for our current considerations. 
Diagnoses (diseases) are modeled as a priori independent binary variables causing a set of 
symptoms (findings) which constitute the bottom layer. The Promedas network currently 
consists of approximately 2000 diagnoses and 1000 findings. 

The interaction between diagnoses and findings is modeled with a noisy-OR structure. 
The conditional probability of the finding given the parents is modeled by n + 1 numbers, 
n of which represent the probabilities that the finding is caused by one of the diseases and 
one that the finding is not caused by any of the parents. 

The noisy-OR conditional probability tables with n parents can be naively stored in 
a table of size 2 n . This is problematic for the Promedas networks since findings that are 
affected by more than 30 diseases are not uncommon. We use efficient implementation of 
noisy-OR relations as proposed by (Takinawa and D'Ambrosio, 1999) to reduce the size of 
these tables. The trick is to introduce dummy variables s and to make use of the property 

OR(z|yi, y 2 , 2/3) = ^2 OR(x|yi, s)OR(s\y 2 , 2/3) (13) 

s 

The interaction potentials on the right hand side involve at most three variables instead 
of the initial four (left). Repeated application of this formula reduces all tables to three 
interactions maximally. 

When a patient case is presented to Promedas, a subset of the findings will be clamped 
and the rest will be undamped. If our goal is to compute the marginal probabilities of the 
diagnostic variables only, the undamped findings and the diagnoses that are not related to 
any of the clamped findings can be summed out of the network as a preprocessing step. The 
clamped findings cause an effective interaction between their parents. However, the noisy- 
OR structure is such that when the finding is clamped to a negative value, the effective 
interaction factorizes over its parents. Thus, findings can be clamped to negative values 
without additional computation cost (Jaakkola and Jordan, 1999). 

The complexity of the problem now depends on the set of findings that is given as 
input. The more findings are clamped to a positive value, the larger the remaining network 
of disease variables and the more complex the inference task. Especially in cases where 
findings share more than one common possible diagnosis, and consequently loops occur, 
the model can become complex. We illustrate some of the graphs that result after pruning 
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Figure 10: Examples of graph structures, corresponding to patient cases generated with one 
disease, after removal of undamped findings and irrelevant disease variables and 
the introduction of dummy variables. Left and right graphs corresponds to an 
"easy" and a "difficult" case respectively. 



of undamped findings and irrelevant diseases and the introduction of dummy variables for 
some patient cases in Figure 10. 

We use the Promedas model to generate virtual patient data by first clamping one 
disease variable to a positive value and then clamping a finding to its positive value with 
probability equal to the conditional distribution of the findings given this positive disease. 
The union of all positive findings thus obtained constitute one patient case. For each patient 
case, the corresponding truncated graphical model is generated. Note that the number of 
disease nodes in this graph can be large and hence loops can be present. 

In this subsection, as well as comparing TLSBP with CVM, we also use another loop 
correction approach, loop corrected belief propagation (LCBP) (Mooij and Kappen, 2007), 
which is based on the cavity method and also improves over BP estimates. We use the 
following parameters for TLSBP: S = 100, M = 5, and no bound b. Again, we apply the 
junction tree method to obtain exact marginals and compare the different errors. Figure 
11 shows results for 146 different random instances. 

We first analyze the TLSBP results compared with BP (Figure 11a). The region in 
light gray color indicates TLSBP improvement over BP. The observed results vary strongly 
because of the wide diversity of the particular instances, but we can basically differentiate 
two scenarios. The first set of results include those instances where the BP error is corrected 
almost up to machine precision. These patient cases correspond to graphs where exhaustive 
enumeration is tractable, and TLSBP found almost all the generalized loops. These are 
the dots appearing in the bottom part of Figure 11a, approximately 14% of the patient 
cases. Note that even for errors of the order of 10 -2 the error was completely corrected. 
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Figure 11: Results of 146 random patient cases with one disease, (a) TLSBP error versus 
BP error, (b) CVM error versus BP error, (c) LCBP error versus BP error, 
and (d) LCBP error versus TLSBP error. 



Apart from these results, we observe another group of instances where the BP error was not 
completely corrected. These cases correspond to the upper dots of Figure 11a. The results 
in these patient cases vary from no significant improvements to improvements of four orders 
of magnitude. 

Figure lib shows the performance of CVM considering all maximal factors together 
with all loops that consist up to four different variables as outer regions. We can see that, 
contrary to TLSBP, CVM in this domain performs poorly. For only one instance the CVM 
result is significantly better than BP. Moreover, the computation time required by CVM 
was much larger than TLSBP in all instances (data not shown). These results can be 
complemented with the study developed in (Mooij and Kappen, 2007), where it is shown 
that CVM does not perform significantly better for other choices of regions. 
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Figure 12: Results of applying TLSBP to 146 patient cases with one disease, (a) Relation 
between computational cost needed by LCBP and TLSBP. (b) Ratio between 
TLSBP and BP errors versus number of loops found, (c) Ratio between TLSBP 
and BP errors versus computation time. 



Figure 11c shows results of LCBP, the approach presented in (Mooij and Kappen, 2007) 
on the same set of instances. As in the case of TLSBP, LCBP significantly improves over BP. 
A comparison between both approaches is illustrated in Figure lid, where those instances 
where TLSBP is better are marked in light gray color. For 41% of the cases TLSBP improves 
the LCBP results, sometimes notably. TLSBP enhacements were made at the cost of more 
time, as Figure 12a suggest, where in 85% of the instances TLSBP needs more time. 

To analyze the TLSBP results in more depth we plot the ratio between the error obtained 
by TLSBP and the BP error versus the number of generalized loops found and the CPU time. 
From Figure 12b we can deduce that cases where the BP error was most improved, often 
correspond to graphs with a small number of generalized loops found, whereas instances 
with highest number of loops considered have minor improvements. This is explained by 
the fact that some instances which contained a few loops were easy to solve and thus the 
BP error was significantly reduced. An example of one of those instances corresponds to 
the Figure 10 (left). On the contrary, there exist very loopy instances where computing 
some terms was not useful, even if a large number of them (more than one million) where 
considered. A typical instance of this type is shown in Figure 10 (right). The same argument 
is suggested by Figure 12c where CPU time is shown, which is often proportional to the 
number of loops found. 

In general, we can conclude that although the BP error was corrected in most of the 
instances, there were some cases in which TLSBP did not give significant improvements. 
Considering all patient cases, the BP error was corrected in more than one order of magni- 
tude for more than 30% of the cases. 
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6. Discussion 

We have presented TLSBP, an algorithm to compute corrections to the BP solution based on 
the loop series expansion proposed by Chertkov and Chernyak. 2 In general, for cases where 
all loops can be enumerated the method computes the exact solution efficiently In contrast, 
if exhaustive enumeration is not tractable, the BP error can be reduced significantly. The 
performance of the algorithm does not depend directly on the size of the problem, but on 
how loopy the original graph is, although for larger instances it is more likely that more 
loops are present. 

We have also shown that the performance of TLSBP strongly depends on the degree 
of coupling between the variables. For weak couplings, errors in BP are most prominently 
caused by small simple loops and truncating the series is very useful, whereas for strongly 
coupled variables, loop terms of different sizes tend to be of the same order of magnitude, 
and truncating the series can be useless. For those difficult cases, BP convergence is also 
difficult, and magnetizations at the fixed point tend to be close to extremum values, causing 
numerical difficulties in the calculation of the loop expansion formulas. In general, we can 
conclude that the proposed approach is useful in an intermediate regime, where BP results 
are not very accurate, but BP is still converging. 

We have confirmed empirically that there is a correlation between the BP result and the 
potential improvements using TLSBP. Those cases where the BP estimate is most corrected 
correspond often to cases where the BP estimate is already accurate. Whether a given BP 
error is acceptable or not depends on the inference task and the specific domain. 

The proposed approach has been compared with CVM selecting loops of four variables 
and maximal factors as outer clusters. For highly regular domains with translation invari- 
ances such as square grids, CVM performs better than TLSBP in difficult instances (strong 
interactions). This is not surprising, since CVM exploits the symmetries on the original 
graph. However, for other domains such as random graphs, or medical diagnosis, TLSBP 
show comparable, or even better results than CVM with our choice of clusters. 

The TLSBP algorithm searches the graph without considering information accessible 
from the BP solution, which is used to compute the loop corrections only as a final step. 
Therefore, it can be regarded as a blind search procedure. We have also experimented 
with a more "heuristic" algorithm where the search is guided in some principled way. Two 
modifications of the algorithm have been done in that direction: 

1. One approach consisted in modifying the third step in a way that, instead of applying 
blind mergings, generalized loops which have larger contributions (largest |r(C)|) were 
merged preferentially. In practice, this approach tended to check all combinations of 
small loops which produced the same generalized loop, causing many redundant merg- 
ings. Moreover, the cost of maintaining sorted the "best" generalized loops caused a 
significant increase in the computational complexity. This approach did not produce 
more accurate results neither was a more efficient approach. 

2. Also, instead of pruning the DFS search for complex loops using the parameter M, we 
have used the following strategy: we computed iteratively the partial term of the loop 

2. The source code of the algorithm and a subset of the datasets used in the experimental section can be 
downloaded from : http://www.cns.upf.edu/vicent/. 
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that is being searched, such that at each DFS step one new term using Equations (6) 
and (7) is multiplied with the current partial term. If at some point, the partial term 
was smaller than a certain threshold A, the DFS was pruned. This new parameter 
A was then used instead of M and result in an appropriate strategy for graphs with 
weak interactions. For cases where many terms of the same order existed, a small 
change of A caused very different execution times, and often too deep searches. We 
concluded that using parameter M is a more suitable choice in general. 

TLSBP can be easily extended in other ways. For instance, as an anytime algorithm. 
In this context, the partition sum or marginals can be computed incrementally as more 
generalized loops are being produced. This allows to stop the algorithm at any step and 
presumably, the more time used, the better the solution. The "improvement if allowed more 
time" can be a desirable property for applications in approximate reasoning, (Zilberstein 
and Russell, 1996). Another way to extend the approach is to consider the search for loops 
as a compilation stage of the inference task, which can be done offline. Once all loops are 
retrieved and stored, the inference task would require much less computational cost to be 
performed. 

During the development of this work another way of selecting generalized loops has been 
proposed (Chertkov and Chernyak, 2006) in the context of Low Density Parity Check codes. 
Their approach tries to find only a few critical simple loops, related with dangerous noise 
configurations that lead to Linear Programing decoding failure, and use them to modify the 
standard BP equations. Their method shows promising results for the LDPC domain, and 
can be applied to any general graphical model as well, so it would be interesting to compare 
both approaches. 

There exists another type of loop correction methods that improves BP, which is quite 
different from the approach discussed here (Montanari and Rizzo, 2005; Parisi and Slanina, 
2006; Mooij et al., 2007; Mooij and Kappen, 2007). Their argument is based on the cavity 
method. BP assumes that in the absence of variable i, the distribution of its Markov blanket 
factorizes over the individual variables. In fact, this assumption is only approximately true, 
due to the loops in the graph. The first loop correction is obtained by considering the 
network with variable i removed and estimating the correlations in the Markov blanket. This 
argument can be applied recursively, yielding the higher order loop corrections. Whereas 
TLSBP computes exactly the corrections of a limited number of loops, the cavity based 
approach computes approximately the corrections due to all loops. An in-depth comparison 
of the efficiency and accuracy of these approaches should be made. 

As a final remark, we mention the relation of the loop series expansion with a similar 
method originated in statistical physics, namely, the high-temperature expansion for Ising 
models. This expansion of the partition function is similar to the one proposed by Chertkov 
and Chernyak, in the sense that every term has also a direct diagrammatic representation 
on the graph, although not in terms of generalized loops. Note however, that the loop ex- 
pansion is relative to the BP result, contrary to the high-temperature expansion. Another 
difference is that the high temperature expansion is an expansion in a small parameter (the 
inverse temperature), whereas the loop expansion has no such small parameter. Finally, an- 
other related approach is the walk-sum framework for inference in certain Gaussian Markov 
Models (Malioutov et al., 2006), where means and covariances between any two nodes of 
the graph have an interpretation in terms of an expansion of walks in the graph. They also 
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show that Gaussian loopy BP has a walk-sum interpretation, computing all walks for the 
means but only a subset of walks for the variances. 

Acknowledgments 

We would like to thank the reviewers for their constructive suggestions that helped us to 
improve the paper. We also acknowledge financial support from the Catedra Telefonica 
Multimedia, the Interactive Collaborative Information Systems (ICIS) project (supported 
by the Dutch Ministry of Economic Affairs, grant BSIK03024), and the Dutch Technology 
Foundation (STW). Finally, we thank Bastian Wemmenhove for providing patient cases of 
the Promedas medical system, and Andreas Kaltennbruner for useful suggestions. 

References 

M. Chertkov and V. Y. Chernyak. Loop series for discrete statistical models on graphs. 
Journal of Statistical Mechanics: Theory and Experiment, 2006(06) :P06009. URL 
http : //arxiv . org/abs/cond-mat/0601487. 

M. Chertkov and V. Y. Chernyak. Loop calculus helps to improve belief propagation and 
linear programming decodings of LDPC codes. In invited talk at 44th Allerton Conference, 
September 2006. URL http://www.arxiv.org/abs/cs/0609154. 

G. Elidan, I. McGraw, and D. Koller. Residual belief propagation: Informed scheduling 
for asynchronous message passing. In Proceedings of the Twenty-second Conference on 
Uncertainty in AI (UAI), Boston, Massachussetts, July 2006. 

W. T. Freeman, E. C. Pasztor, and O. T. Carmichael. Learning low-level vision. Int. J. 
Comp. Vision, 40:25-47, 2000. 

R. G. Gallager. Low-density parity check codes. MIT Press, 1963. 

T. Heskes, K. Albers, and H. J. Kappen. Approximate inference and constraint optimisation. 
In Proceedings UAI, pages 313-320, 2003. 

T. Jaakkola and M. I. Jordan. Variational probabilistic inference and the QMR-DT network. 
Journal of artificial intelligence research, 10:291-322, 1999. 

D. B. Johnson. Finding all the elementary circuits of a directed graph. SI AM J. Comput., 
4(l):77-84, 1975. 

M. Jordan, Z. Ghahramani, T. S. Jaakkola, and L. Saul. An introduction to variational 
methods for graphical models. In Learning in Graphical Models, pages 105-161. Cam- 
bridge, MA: MIT Press, 1999. 

F. R. Kschischang, B. J. Frey, and H.-A. Loeliger. Factor graphs and the sum-product 
algorithm. IEEETIT: IEEE Transactions on Information Theory, 47, 2001. 



29 



ViCENg Gomez and J.M. Mooij and H.J. Kappen 



S. L. Lauritzen and D. J. Spiegelhalter. Local computations with probabilities on graphical 
structures and their application to expert systems. J. Royal Statistical society B, 50: 
154-227, 1988. 

M. A. R. Leisink and H. J. Kappen. A tighter bound for graphical models. Neural Comput., 
13(9):2149-2171, 2001. 

D. Malioutov, J. Johnson, and A. Willsky. Walk-sums and belief propagation in gaussian 
graphical models. Journal of Machine Learning Research, 7, October 2006. 

R. McEliece, D. MacKay, and J. Cheng. Turbo decoding as an instance of Pearl's belief 
propagation algorithm. Journal of Selected Areas of Communication, 16:140-152, 1998. 

M. Mezard, G. Parisi, and R. Zecchina. Analytic and algorithmic solution of random 
satisfiability problems. Science, 297:812-815, 2002. 

A. Montanari and T. Rizzo. How to compute loop corrections to the Bethe approximation. 
Journal of Statistical Mechanics: Theory and Experiment, 2005(10):P10011, 2005. URL 
http : //arxiv . org/abs/cond-mat/0506769. 

J. M. Mooij and H. J. Kappen. Loop corrections for approximate infer- 
ence. Journal of Machine Learning Research, 8:1113-1143, May 2007. URL 
http : //arxiv . org/abs/cs/0612030. 

J. M. Mooij and H. J. Kappen. On the properties of the Bethe approximation and loopy 
belief propagation on binary networks. Journal of Statistical Mechanics: Theory and 
Experiment, 2005(11):P11012, 2005. 

J. M. Mooij, B. Wemmenhove, H. J. Kappen, and T. Rizzo. Loop corrected belief propa- 
gation. In Proceedings of the Eleventh International Conference on Artificial Intelligence 
and Statistics, AISTATS, 2007. 

K. P. Murphy, Y. Weiss, and M. I. Jordan. Loopy belief propagation for approximate 
inference: An empirical study. In Proceedings of Uncertainty in AI, pages 467-475, 1999. 

J. Parisi and F. Slanina. Loop expansion around the Bethe-Peierls approximation for lattice 
models. Journal of Statistical Mechanics, page L02003, 2006. 

J. Pearl. Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference. 
Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1988. ISBN 1558604790. 

A. Pelizzola. Cluster variation method in statistical physics and probabilistic graphical 
models. Journal of Physics A: Mathematical and General, 38(33):R309-R339, 2005. URL 
http : //arxiv . org/abs/cond-mat/0508216. 

G. Potamianos and J. Goutsias. Stochastic approximation algorithms for partition function 
estimation of Gibbs random fields. Information Theory, IEEE Transactions on,, 43, 6: 
1948-1965, Nov 1997. 



30 



Truncating the Loop Series Expansion for BP 



M. A. Shwe, B. Middleton, D. E. Heckerman, M. Henrion, E. J. Horvitz, H. P. Lehman, 
and G. F. Cooper. Probabilistic diagnosis using a reformulation of the internist-1/ QMR 
knowledge base. Methods of Information in Medicine, 30:241-55, 1991. 

J. Sun, Y. Li, S. B. Kang, and H. Y. Shum. Symmetric stereo matching for occlusion 
handling. Proceedings CVPR, 2:399-406, 2005. 

M. Takinawa and B. D'Ambrosio. Multiplicative factorization of noisy-MAX. In Proceedings 
of the 15th Conference on Uncertainty in Artificial Intelligence, pages 622-30, 1999. 

R. E. Tarjan. Enumeration of the elementary circuits of a directed graph. SI AM J. Comput., 
2(3):211-216, 1973. 

J. C. Tiernan. An efficient search algorithm to find the elementary circuits of a graph. 
Commun. ACM, 13(12):722-726, 1970. 

M. Wainwright, T. Jaakkola, and A. Willsky. A new class of upper bounds on the log 
partition function. 51 (7):2313-2335, July 2005. 

M. Welling, T. Minka, and Y. W. Teh. Structured region graphs: Morphing EP into GBP. 
In Proceedings of the 21th Annual Conference on Uncertainty in Artificial Intelligence 
(UAI-05), page 609, Arlington, Virginia, 2005. AUAI Press. 

W. Wiegerinck, H. J. Kappen, E. W. M. T. ter Braak, W. J. P. P. Burg, M. J. Nijman, 
Y. L. O, and J. P. Neijt. Approximate inference for medical diagnosis. Pattern Recognition 
Letters, 20:1231-1239, 1999. 

J. S. Yedidia, W. T. Freeman, and Y. Weiss. Generalized belief propagation. In T.K. 
Leen, T.G. Dietterich, and V. Tresp, editors, Advances in Neural Information Processing 
Systems 13 (Proceedings of the 2000 Conference), 2001. 

J. S. Yedidia, W. T. Freeman, Y. Weiss, and A. L. Yuille. Constructing free-energy ap- 
proximations and generalized belief propagation algorithms. Information Theory, IEEE 
Transactions on, 51:2282-2312, 2005. 

A. L. Yuille. CCCP algorithms to minimize the Bethe and Kikuchi free energies: Convergent 
alternatives to belief propagation. Neural computation, 14:1691-1722, 2002. 

S. Zilberstein and S. Russell. Approximate reasoning using anytime algorithms. AI maga- 
zine, 17(3):73-83 (lp.1/4), 1996. ISSN 0738-4602. 



31 



